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SUMMARY 


A  computational  fluid  dynamics  code  named  MULTIFLOW-2D  has  been  acquired  and  is 
being  put  into  service  for  the  numerical  simulation  of  gas  turbine-type  combustors. 

Certain  areas  of  the  code  were  found  to  be  in  need  of  improvement  for  this  application; 
these  are: 

•  The  grid  generation  algorithm 

•  The  manner  of  discretisation  of  the  continuity  and  momentum  equations,  and 

•  The  treatment  of  the  inlet  boundary  conditions. 

Improvements  to  these  areas  were  developed  and  incorporated  to  the  code,  and  the  code 
functioning  in  non-reactive  flows  was  tested  by  means  of  numerical  experiments;  verification 
against  physical  experiments  is  planned. 
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1.  INTRODUCTION 


As  part  of  a  program  to  develop  a  capability  in  computational  modeling  of  gas  turbine  com¬ 
bustor  flows,  the  MULTIFLOW-2D  code  is  being  adapted  for  application  to  this  purpose. 
The  code  was  obtained  from  the  Aero  Propulsion  and  Power  Directorate  of  the  US  Air 
Force  Wright  Laboratory  under  the  Mutual  Weapons  Development  Data  Exchange  Agree¬ 
ment  Annex  7017. 

The  code  has  been  modified  both  to  correct  deficiencies  and  to  introduce  features  necessary 
for  the  intended  application. 

2.  THE  MULTIFLOW.2D  CODE 

2.1  Brief  description 

The  theoretical  basis  of  this  code  is  detailed  by  its  author  in  Refs  [1]  and  [2],  while  further 
description  of  the  code  capabilities  and  general  features  may  be  found  in  [3]. 

In  brief,  the  code  finds  solutions  to  the  conservation  equations  of  flow  discretised  in  finite 
volume  form.  The  discretisation  points  form  a  grid  which  covers  the  flow  region  using  an 
equal  number  of  points  per  grid  line  (Boundary  Fitted  Coordinate  system).  The  flow  may 
be  two-dimensional  or  axisymmetric  (including  axisymmetric  rotation,  termed  swirl). 

Solutions  are  found  for  the  conservation  of  mass,  linear  momentum,  energy  and  chemical 
species.  Closure  of  the  conservation  equations  for  turbulent  flows  is  obtained  by  means 
of  the  two-equation  k  —  t  model;  the  stochastic  turbulent  diffusion  flame  model  uses  an 
assumed  probability  distribution  function  generated  using  its  first  two  moments  (mean  and 
rms),  which  are  obtained  from  their  own  transport  equations.  There  is  also  a  facility  to 
simulate  liquid  fuel  sprays. 

The  solution  process  includes  a  multigrid  system,  whereby  coarse  grids  are  utilised  to  remove 
the  low  frequency  errors  which  slow  down  the  convergence  of  the  fine  grids. 

2.2  Improvements  required 

The  code  was  installed  in  the  local  computer  system,  “debugged”  iuid  tested  until  qualita¬ 
tively  satisfactory  computations  were  obtuned  for  a  number  of  cases  including  that  of  an 
annular  jet  discharging  into  a  sudden  expansion  (with  and  without  a  coflowing  centre  jet) 
and  a  separating  diffuser.  However,  the  operation  of  certain  code  functions  was  deemed 
unsatisfactory,  namely; 


•  Great  difficulty/inability  to  use  fine  meshes,  and 

•  Great  difficulty /inability  to  compute  flows  with  variable  density. 


Having  achieved  a  degree  of  confidence  in  the  operation  of  the  code  development  work  was 
initiated  to  overcome  the  observed  difficulties.  Section  3  of  this  Report  describes  the  addi¬ 
tions  emd  modifications  introduced  to  this  effect.  Difficulties  in  using  very  fine  meshes  were 
attributed  to  the  manner  in  which  coarse  grids  were  generated  from  the  finest  grid,  which 
is  discussed  in  Section  3.1,  and  to  the  manner  in  which  the  coefficients  of  the  discretised 


equations  were  generated,  which  is  discussed  in  Section  3.3.2.  Difficulties  in  the  compu¬ 
tation  of  flows  with  variable  density  were  attributed  to  the  need  for  a  more  sophisticated 
discretisation  of  the  continuity  equation  to  take  into  account  aero-thermodynamic  relations, 
which  is  discussed  in  Section  3.2,  and  to  the  need  to  reformulate  inlet  boundary  conditions 
to  conserve  total  pressure,  which  is  discussed  in  Section  3.3.1. 

Validation  of  the  new  code  against  experimental  results  in  non-reacting  flows  is  planned  for 
the  near  future.  Section  4  presents  examples  of  computational  outputs  to  demonstrate  the 
functioning  of  the  code  and  some  noteworthy  features. 

3.  DEVELOPMENT  OF  MULTIFLOW-2D 
3.1  Grid  generation  process 

The  difficulties  experienced  using  very  fine  grids  were  investigated,  and  were  attributed  to 
the  practice  of  generating  fine  grids  by  subdivision  of  coarse  grids.  This  results  in  ill-fitting 
meshes,  particularly  about  finite-radius  comers,  and  is  unsuitable  for  continuous  curved 
boundaries,  as  shown  in  Fig  1.  The  grid  generation  algorithm  was  therefore  turned  about, 
to  fit  the  finest  mesh  to  the  boundaries  and  obtain  coarse  meshes  by  skipping  every  second 
grid  line,  as  illustrated  in  Fig  2.  This  allowed  some  examples  to  be  computed  to  maximum 
grid  fineness  (6  grid  levels),  which  was  previously  unachievable. 

The  revised  scheme  may  introduce  error  in  the  process  of  changing  from  a  coarse  to  a  fine 
grid  and  vice  versa.  However,  no  error  is  introduced  in  the  transformation  of  velocities  on 
the  basis  of  conservation  of  mass.  Fig  3  shows  that  the  projected  areas  dA  are  not  affected 
by  the  fact  that  a  coarse  cell  does  not  fully  enclose  four  fine  cells,  because  the  extreme 
corners  are  common.  It  is  still  exact  to  say  that,  for  instance: 

(pUdA)^  =  (pUdA)fj  ipUdA)fj.,  (1) 

where  p  is  density  and  U  the  component  of  velocity  normal  to  A. 

On  the  other  hand,  scalar  quantities  are  transformed  on  the  assumption  that  a  fine  cell  is 
1  /4  the  size  of  a  coarse  cell,  and  thence  the  centre  of  the  coarse  cell  is  the  CG  of  the  centres 
of  the  fine  cells.  Fig  3  shows  that  in  the  revised  scheme  an  error  may  be  introduced  in  that 
the  point  at  which  the  coarse  value  is  computed  is  C,  where  the  coarse  value  of  a  scalar 
quantity  5  is  given  by: 


sfj = lisfj  -b  sr_,j  -b  sfj., + (2) 

However,  the  coarse  quiuitity  5,^  should  have  been  computed  not  at  C  but  at  the  coarse 
cell  centre  Cfj. 

This  added  inaccuracy  in  the  transfer  of  scalars  between  grids  is  of  no  consequence  since 
in  the  Multigrid  system  used  in  MULTIFL0W-2D  the  finest  grid  (the  reported  solution 
grid)  must  in  the  end  be  converged  on  its  own.  Simple  and  relatively  inaccurate  transfer 
mechanisms  between  grids  have  no  effect  on  the  end  result  on  the  finest  grid. 
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*; 


•  •  • 


•  • 


3.2  Discretisation  of  the  Continuity  equation 


In  the  original  formulation  of  MULTIFLOW-2D  the  continuity  equation  was  discretised 
as  for  constant-density  Oow.  For  use  in  combustion  it  is  necessary  to  be  able  to  compute 
flows  with  large  density  variations,  and  where  density  is  linked  to  temperature  and  pressure 
according  to  the  laws  of  aero-thermodynamics. 

Two  areas  require  development  in  this  respect:  evaluation  of  the  density-velocity  product 
at  a  cell  face,  and  compliance  with  the  second  law  of  thermodynamics. 

3.2.1  Density  interpolation 

The  mass  flow  rate  across  a  cell  face  dA  may  be  expressed  as: 


pUdA.  (3) 

In  the  staggered-grid  system  presently  used  all  scalar  quantities  (e.g.,  pressure  and  density) 
are  computed  at  the  cell  centre.  A  cell  is  identified  by  its  t,ji  indices  and  the  coordinates 
of  the  top  right-hand  corner.  Vectorial  quantities  U  and  V  are  computed  at  the  centre  of 
the  cell  faces  adjacent  to  the  identifying  comer  (Fig  4).  It  thence  becomes  necessary  to 
interpolate  density  to  obtain  its  value  at  the  two  cell  faces  for  U  and  V.  The  challenge  is 
to  And  an  interpolation  algorithm  which  is  computationally  inexpensive  and  has  a  basis  in 
the  physics  of  the  problem. 

One  possible  solution  is  standard  spline  fitting,  in  itself  a  computationally  expensive  pro¬ 
cess.  Moreover,  this  requires  four  data  points  straddling  each  cell  face,  which  leads  to  a 
computational  stencil  of  at  least  9  cells  (Fig  5).  This  gives  rise  to  a  number  of  special  cases 
near  boundaries,  a  situation  to  be  avoided  as  it  is  also  computationally  expensive. 

At  the  other  extreme  the  simple  arithmetic  average  between  adjacent  cell  centres  may  be 
adopted.  This  is  computationally  inexpensive  and  only  requires  the  minimal  5-cell  stencil. 
On  the  other  hand,  no  information  is  transmitted  regarding  the  physics  of  the  problem,  i.e., 
upstream  and  downstream  values  and  gradients  have  the  same  weight. 

An  alternative  has  been  given  by  Karki  and  Patankar  [4],  who  adopt  an  upwinding  approach, 
illustrated  in  Fig  6  for  the  west  face.  Calling 


M{<f>)  =  Largest  of  <f>  and  0  , 
then  their  interpolated  value  is: 


M(U„) 

'’-'"“iAr 


+  Pp 


t/u, 


(4) 


(5) 


As  a  result,  density  on  the  west  face  will  be  the  density  at  the  west  cell  centre  pw  if  the 
flow  is  conventionally  positive  (West  to  East),  and  will  be  the  density  at  the  cell  centre 


Pp  if  the  flow  is  from  East  to  West.  The  underlying  physical  assumption  is  that  density  is 
primarily  a  convected  property. 

When  coarse  cells  are  utilised  luad  combustion  is  expected  it  is  desirable  to  retain  some 
links  to  adjacent  values.  To  this  end  the  upwind  concept  has  been  extended  to  achieve  a 
form  which  is  as  computationally  simple  as  the  last  two  above  but  still  carries  some  higher 
level  information.  This  is  achieved  by  firstly  considering  an  upwind-biased  second  order 
polynomial  fit,  i.e.,  with  two  data  points  upwind  of  the  cell  face  and  one  downwind.  This 
is  exemplified  in  Fig  7  for  the  computation  of  density  at  the  p  face  (the  face  corresponding 
with  Up). 

Elementary  algebra  yields  for  a  conventionally  positive  U : 

Pp  = -^(^Pe  +  6pp  -  pw)  ■  (6) 

For  U  of  the  opposite  sense: 

Pp  =  +  ^Pe  -  pee)  ■  (7) 

The  switching  function  M(U)  may  be  used  now  to  select  one  of  the  interpolating  functions. 
With  previous  notation  the  interpolated  values  at  the  four  faces,  as  used  in  Fig  8,  are: 


Pp  -  [(3p£:  +  6/}p  -  pw)M{Up)  4-  (3/»p  +  6pE  -  PEE)M{-Up)]  (at  Up)  ,  (8) 

Pw  =  ((3pp  +  6pw  —  Pww)M(Uu,)  +  (3^^  -4  6pp  —  Pb)M{-Uw)]  (at  U^,) ,  (9) 

Pp  =  ^  l(3p^  -I-  6pp  -  ps)M(Vp)  -t-  (3pp  -I-  6pN  -  Pnn)M(-Vp)]  (at  Vj,) ,  (10) 

p,  =  ^  [(3pp  -H  6ps  -  pss)M{V,)  +  {Zps  6pp  -  Ps)M(-V,)]  (at  V,) .  (11) 

This  is  not  as  computationally  expensive  as  the  third-order  spline  fit  but  again  it  involves 
a  9-cell  stencil  (note  cells  EE,  SS,  NN,  and  WW). 

A  simplified  algorithm  is  adopted  instead,  which  may  simply  be  described  as: 


Pjace  —  ^{.Pupvnnded  4*  Paveraged)  • 


For  Pp  at  the  Up  face  and  West  to  East  flow  the  above  expression  yields: 


(12) 


1/-  ,  Pp  +  Pe^ 

Pp  -  -^{PP  +  — ^)  • 


(13) 


On  the  other  hand  the  upwind-biased  parabolic  interpolation  gives: 


(14) 


Pp  -  g(3p£:  +  6pp  —  pw)  • 


The  difference  is: 


^Pp  ~  -^{PE  -  Pw)  ■ 


(15) 


This  is  almost  an  order  of  magnitude  smaller  than  the  differtnce  in  densities  between  up¬ 
stream  and  downstream,  i.e.,  the  simplified  algorithm  (12)  is  a  very  close  approximation  to 
the  upwind-biased  parabolic  fit. 

It  is  seen  that  the  simplified  interpolation  algorithm  has  the  following  properties: 


•  Computationally  inexpensive  (involves  only  the  minimal  stencil), 

•  Incorporates  some  upwind  biasing  in  response  to  the  physics  of  the  problem,  and 

•  Incorporates  some  higher  level  information,  i.e.,  upwind  slope. 

3.2.2  Complying  with  the  Second  Law  of  Thermodynamics 


Consider  the  computational  cell  P  and  adjacent  locations  shown  in  Fig  8. 
The  continuity  equation  may  be  written: 


pU  DU  Ip  ~pU  DU  U  -pV  DV  I.  -fpV  DV  |p=  0  ,  (16) 

where  DU,DV  the  cell  face  areas  normal  to  U,V.  Densities  are  assumed  interpolated 
as  per  the  previous  discussion. 

The  exact  values  of  the  variables  are  employed  above,  since  the  sum  is  equated  to  zero. 
Exact  values  can  be  related  to  current,  approximate  values  (indicated  by  an  asterisk)  by 
means  of  corrections  (indicated  by  primes); 


•  • 


•  • 


•  • 


•  •  • 


•  • 


•  • 


p  =  p^  +  p  , 

(17) 

• 

• 

U  =  U’  +  U', 

(18) 

V  =  V*  V'  . 

(19) 

Then,  neglecting  second-order  products: 

• 

• 

pU  =  p'l/'-f  p*l/' 4- p*l/*  , 

(20) 

pV  =  p'V  +  p^V'  +  p^V  . 

(21) 

• 

• 
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•  • 


Substituting  and  rearranging: 


T 


■*; 


p'U'  DU  Ip  -p’V  DU  U  ~p'V‘  DV  |.  +p'V'  DV  |p 
+pU'  DU  Ip  -pU’  DU  U  -pV'  DV  U  +p' V  DV  |p  = 

-  [p'U'  DU  Ip  -p'U'  DU  U  -p'V  DV  I,  +p’V  DV  |p]  .  (22) 


The  right-hand  side  is  the  residual  obtained  when  computing  the  continuity  equation  with 
all  current  values,  which  will  be  called  RC. 


The  first  four  terms  on  the  left-hand  side  involve  current  values  of  density  and  correction  val¬ 
ues  for  velocity.  These  are  kept  as  they  are  and  will  be  used  to  compute  the  four  corrections 


U  U  V'  V' 

'pi  • 


The  next  four  terms  are  a  challenge  because  they  must  be  put  in  a  form  compatible  with 
the  discretised  Navier-Stokes  equations,  i.e.,  in  terms  of  velocity  and  pressure.  The  current 
values  of  velocity  are  left  as  they  are,  and  the  task  is  to  express  the  density  corrections  p  , 
which  are  placed  at  the  f/,V  face  locations,  in  terms  of  the  pressure  correction  pp,  which  is 
located  at  the  cell  centre.  This  must  be  done  in  such  a  manner  that  the  behaviour  of  the 
discretised  equations  is  correct  over  the  range  of  Mach  numbers,  changing  from  elliptic  to 
hyperbolic  as  the  flow  changes  from  subsonic  to  supersonic. 


The  final  form  of  the  discretised  continuity  equation  will  be: 


# 


4r; 


•  • 


•  • 


•  •  • 


*^iU^  +  ujf/p  +  azV,  +  a^Vp  +  OsPp  =  F  , 


(23) 


1 


I 


where: 


•  • 


Oi  =  -p’DU  !«, , 

(24) 

a,  =  p-Z)f/|p, 

(25) 

• 

• 

03  =  -p-DVl, 

(26) 

a,  =  p-nvip  , 

(27) 

and  expressions  must  be  found  for  as  and  F. 

• 

• 

The  first  step  in  linking  p'  and  p  is  to  transform  the  four  terms  into  expressions  containing 
the  cell-centred  values  by  means  of  the  interpolation  formulae.  The  remaining  four  terms, 
i.e.,  the  last  four  terms  of  Eq  (22)  were: 

• 

• 

pU‘  DU  Ip  -pU‘  DU  U  -p'V  DV  |.  -Vp'V'  DV  |p  . 

(28) 

Into  this  expression  substitute  the  interpolated  values  of  p  obtained  using  the  interpolation 
procedure,  i.e.: 

6 

• 

• 

•  •  • 


(29) 


1  ,  M(U,)  1  ,  M(-U,)  + 

2PP—  +  r^—UT'  ’ 

1  /  M([/„)  1  >  M(—UJ)  lp/>  +  Pw,_^  ,,  ^ 

i'^-uT  r'—UT'  ■"  2—2—^“  ■ 

1 .  urn  1 ,  M(-n)  ip,+,;„  . 

2"'—  +  2'’"—^  +  >*‘  ’ 

- y -  +  ^Pp - y -  +  2 - 2 - 


(Asterisks  have  been  conveniently  dropped  since  only  current  values  of  velocity  are 
involved) 

Substituting  Eqns  (29)-(32)  into  (28)  and  rearranging; 


Ppmu,)  +  +  PsiMi-U,)  +  , 

-Pw(M(u^) + -  p'pm-u.) + y)^  ’ 
+/>;.( A/(Vp)  +  ^)^  +  Pr,(M(-V,)  +  ^)^  , 
~Ps(A^(Vi)  +  ^^“2 - +  '^)“2~  ’ 

Define  and  collect  coefficients  VA  (for  ‘Velocity  times  Area’): 


Terms  for  pg  :VA£  =  (M(-Up)  +  y)^^  . 

V  DV 

Terms  for  p'^  :VAn  =  {M(-V^)  +  , 

V  DV 

Terms  for  p's  :VAs  =  -(A/(  V.)  +  ^)-^  , 
Terms  for  p^  :VA»v  =  -(M(  U,)  +  ^)^  , 


and  for  pp: 


VAp  =  (M(U,)+^)^-(M{-U^)+^)^+{M{V,)A^)^-{M{-V,)+^)^  . 


The  four  terms  (28)  can  now  be  expressed  as  five  terms  in  cell-centred  density  corrections: 


7 


V Appp  +  V AePe  +  ^ AsPs  +  ^ AsPs  +  ^ Awpw  ■  (39) 


The  next  step  is  to  relate  the  cell-centred  density  corrections  to  cell-centred  pressure  cor¬ 
rections.  From  the  state  equation: 


dp  dp  dt 

P  P  t  ' 


(40) 


and  from  the  definition  of  entropy: 


dt  da  —  I  dp 

t  -y  P  ' 


where  s  is  entropy  per  unit  mass  in  units  of  Cp  is  the  specific  heat  at  constant 

pressure  and  7  the  adiabatic  exponent. 


•  • 


•  • 


•  • 


Eliminating  dt/t  and  using  the  definition  of  isentropic  velocity  of  sound 


c*  =  7P/P  , 


(42) 


•  •  • 


results  in: 


dp=^(dp-^da)  .  (43) 

This  equation  can  be  put  in  Eulerian  form  as: 

Dp  =  ^{Dp-^Da).  (44) 

Substituting  the  corrections  for  the  differentials: 

p'  =  j,{p-^^{Da)'),  (45) 

where  (Da)'  is  a  correction  to  the  process  entropy  increase,  arising  from  the  corrections  to 
velocities,  densities,  etc. 


•  • 


•  • 


•  • 
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•  •  • 


Call  VH'{{ot  ‘Viscous  Heating’): 


(46) 


and  Eq  (45)  is  written: 


p  =  i(p'  -  VH‘)  .  (47) 

In  isentropic  flow  a  pressure  change  p'  produces  a  density  change  p'  =  p  /c^.  Viscosity 
raises  the  gas  temperature,  diminishing  a  density  increase  (if  p'  >  0)  or  reinforcing  a  density 
reduction  (if  p  <  0).  Note  that  it  is  always  VH'  >  0. 

Substituting  p'  the  five  terms  (39)  can  be  manipulated  to  yield: 


if+  £  (^p)-  E 

E,N,S,W  P,E,N,S,W  ^ 


This  yields  the  desired  coefficient  of  pp  as: 


(48) 


VA 


(49) 


Terms  in  the  surrounding  pressure  corrections  ^p'  are  moved  to  the  right-hand  side,  to  be 
substracted  from  the  continuity  residual;  thus,  so  far: 


F  =  RC-  (— p').  (50) 

E,N,S,W  ^ 

At  this  stage  it  is  possible  to  consider  the  behaviour  of  the  discretised  equation.  Neglecting 
residuals  and  entropy  terms  the  continuity  equation  (from  Eqns  (23)-(27)  and  (49))  is  of 
the  form: 


pU'DU  +  ^^p.  (51) 

This  can  be  rearranged  as: 

^{pUU'  +  M^p).  (52) 

The  relative  importance  of  the  pressure  and  velocity  corrections  is  seen  to  be  proportional 
to  the  square  of  the  Mach  number.  At  low  Mach  numbers  the  velocity  correction  term  will 
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Nw.  '' 


•  • 


•  • 


•  •  • 


•  • 


] 


•  •  • 


predominate  as  is  appropriate  for  elliptic  flows,  whereas  at  high  Mach  numbers  the  pressure 
term  becomes  dominant  as  is  appropriate  for  hyperbolic  flows. 

There  remains  to  express  VH‘  in  terms  of  velocity  and  pressure  corrections.  This  is  done 
by  means  of  the  energy  equation  taken  in  the  following  form  [5]: 


Ds  k 

p—  =  —  +  +  Radiation  .  (53) 

Ut  t 

Here,  k  is  the  coefficient  of  heat  conduction,  r  is  time  and,  with  F  the  effective  viscosity, 
the  viscous  dissipation  function  $  is  given  in  tensorial  form  by  [5]: 


$  =  r 


dui  2  2 


(54) 


where 


X  =  div  u  . 


(55) 


Hence 


$ 

£)s  =  — £)r  +  conduction  +  radiation.  (56) 

pt 

Now,  Da  is  the  entropy  increase  due  to  the  mechanical  work  of  the  viscous  forces,  called 
$.  For  the  present  purposes  what  is  needed  is  the  change  in  the  entropy  increase  due  to 
small  corrections  to  U  and  V  which  produce  small  corrections  in  $  (neglect  the  effect  of 
density  corrections).  In  this  case  the  conduction  and  radiation  terms  will  be  substantially 
unchanged  and  the  entropy  correction  may  be  approximated  by: 


Ws)'  =  ^Dr  . 


Hence,  substituting  (I)^)'  into  VH'  and  manipulating: 


(57) 


VH'  =  (7  -  l)I>r  . 


(58) 


For  Dt,  adopt  the  time  a  parcel  of  fluid  remdns  in  the  reference  volume,  i.e.,  the  transit 
time. 


The  changes  in  $  due  to  U',  V'  are  complex  and  tedious  to  compute.  They  are  obtsuned 
starting  from: 


aVp'^  dv/* 


(59) 
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To  continue  it  is  necessary  to  relate  the  face  velocities  to  the  velocity  components  in  the 
appropriate  coordinate  system.  For  instance,  for  axisymmetric  flow  with  swirl  in  cylindrical 
polar  coordinates  the  expression  for  the  time-averaged  dissipation  function  reduces  to  [6][7]: 


$  =  r[2(t/i  -I-  v;  -t-  {Vlyf)  {U,  -1-  -l-  K  -  Wy)*  -  Ix^]  -  Ipxk  ,  (60) 

where  k  is  the  specific  kinetic  energy  of  the  turbulence,  u>  is  the  swirl  component  of 
velocity  and 


X  =  div  \i  =  Ut  +  Vy  +  V/y  .  (61) 

There  still  remains  to  consider  possible  generalised  coordinates.  This  has  all  been  imple¬ 
mented  but  will  not  be  detailed  further.  Assume  the  derivatives  of  $  required  in  Eq  (59) 
are  known  and  define  coefficients 


ZU  =  ~iT-l)Ardi/dU , 

VA 

ZV  =  -^(n-\)ATd^ldV . 

Hence  the  remaining  viscous  heating  terms  may  be  put  in  the  following  form: 


P.N,S,E,W 


-zvyVy  -  zv.v;  -  zt/pf/;  -  zvjj'y, 

-ZVX  -  ZVpVp  -  ZUnK  -  ZUr^^Uiy, 
-zv.v;  -  zv,.v;,  -  ZU.U',  -  ZU^U'yy, 
-zv.v;  -  ZV„Vi  ~  ZU.U',  -  ZUyU'y 
-ZK,v:  -  zv„,v:  -  ZUy,K.  -  ZUy^U'„, 


This  may  be  efficiently  organised  by  means  of  a  subroutine  which  computes  ZU  and  ZV  at 
a  generic  ceil  location,  and  using  the  subroutine  at  cells  N,  S,  E  and  W  separately.  Then, 
the  computation  at,  for  instance,  N,  would  in  fact  be  the  computation  of  the  second  line  of 
the  above  list;  the  correspondence  is: 


Original:  -ZVX  -  ZV^  “  ZUnU'„  -  ZU„X^  ,  (65) 

Computed  at  N:  -ZV^  -  ZVX  -  ZU^p  -  ZUp,U'„  .  (66) 

That  is,  of  the  values  computed  at  N  keep  —ZV,  (to  be  added  to  the  coefficient  of  V^')  and 
move  all  three  others  to  the  right-hand  side. 
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The  stencil  of  Fig  9  may  assist  in  understanding  this  example.  Note  that  it  is  not  necessary 
to  apply  the  subroutine  at  the  centre  cell. 

The  twelve  terms  ZU  U‘,ZV  V'  added  to  the  RHS  complete  the  desired  function  F.  The 
viscous  heating  components  to  be  added  to  coefficients  Ui  to  are: 


—2ZV,  computed  at  N,  add  to  (i.e.,  for  V^)  ,  (67) 

—2ZUw  computed  at  E,  add  to  aj  (i.e.,  for  U^)  ,  (68) 

—2ZVp  computed  at  S,  add  to  03  (i.e.,  for  V^’)  ,  (69) 

—2ZUt  computed  at  W,  add  to  Oi  (i.e.,  for  U'^)  .  (70) 


Care  must  be  taken  in  regions  close  to  boundaries,  where  U,  V  and  p  are  prescribed  or  •  • 

otherwise  fixed  and  there  must  be  no  computation  of  corrections. 

Recapitulating,  a  routine  is  set  up  to  compute  the  viscous  heating  terms  {ZU,ZV).  This 
routine  is  applied  at  the  four  adjacent  cells  (N,  S,  E,  W)  and  components  found  for  coeffi¬ 
cients  ai  to  and  the  right-hand  side  term  F.  •  • 

Coefficients  oi  to  04  are  completed  by  adding  the  appropriate  density  times  area  components 
(±p*  DU,±p*  DV)  as  per  Eqns.  (24)  to  (27). 

Coefficient  us  is  obtained  as 

•  •  • 


VA 

c3 


\P  . 


(71) 


and  all  surrounding  terms 


•  • 


(72) 


are  moved  to  the  RHS. 

Finally,  the  mass  flow  residual  computed  with  current  values  of  U,  V  and  p  is  added  to  F. 
The  finished  discretised  continuity  equation  is: 

aiU'^  +  +  “3V,'  +  -f  a^pp  =  F  (73) 

This  form  is  suitable  for  all  Mach  numbers  (within  the  limits  of  applicability  of  the  transport 
equations),  and  for  non-isentropic  flows. 


•  • 


•  • 
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•  • 


•  •  • 


3.3  Discretiaing  the  IVaiMport  Equations 

This  Section  considers  improvements  needed  to  the  discretised  transport  equations;  a  brief 
description  of  the  equations  is  given  below. 

The  steady-state  transport  equation  of  a  variable  ^  in  generalised  coordinates  (,Tf,  and 
integrated  o' /er  the  cell  volume  dV  takes  the  form: 


I 

J 


dipU<^)  ^  d(pVd) 


di 


dfj 


dv  =  - 


d  d^  84  a  84  84. 


dV  -t-  S^dV  . 

(74) 


In  these  equations,  U  and  V  are  the  co-variant  transformed  velocities  along  the  transformed 
axes  ^  and  t},  J  is  the  Jacobian  of  the  transformation,  F  the  coefficient  of  diffusion,  S* 
the  source  term  (generation/destruction  of  4  unit  volume  per  unit  time)  and  other 
definitions  are: 


U 

V 

J 

dV 

Cx 

Ca 

Cj 

C, 


8y  8x 

(75) 

8x  8y 

(76) 

8x  8y  8x  8y 

did^i  "  driac 

(77) 

JRAijA( , 

(78) 

^{8xl8rif  {ay! 8rif 

J 

(79) 

^{dxl8if-^(8yl8(.f 
^  J  ' 

(80) 

^{axl8i){8xl8ri)  -I-  {8yl8i){8yl8ri) 

J 

(81) 

C-x. 

(82) 

V2uiable  R  is  unity  in  a  2D  Cartesian  system,  and  equals  the  radial  coordinate  y  in  cylindrical 
polar  coordinates. 


The  source  term  is  usually  assumed  to  be  a  linear  function  of  the  variable: 


S*  =  SU*  +  4SP^ .  (83) 

For  the  transport  of  variables  u  and  v  (i.e.,  velocities  in  the  original  system  of  coordinates 
x,y),  the  equations  become  the  momentum  transport  or  Navier-Stokes  equations  and  the 
ensemble-averaged  source  terms,  expressed  in  x,y  coordinates  for  brevity,  are  (see,  for  in¬ 
stance,  [6][7]): 
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•  • 


•  • 


•  •  • 


•  • 


•  • 


•  • 


•  • 


•  •  • 


5“  = 


dx  dx 


(84) 


(„du\  d  (j,dv\  2d  r  ( du  dv\ 

^  dp  d  fj^du\  d  (r^dv\  2  a  .  2r  f  dv  v\  2p«  piy* 


These  equations  are  discretised  using  finite  differences  techniques  and  applied  at  the  four 
faces  of  a  cell.  This  procedure  results  in  four  Finite  Differences  Equations  (PDE’s)  in  five 
unknowns:  the  four  face  velocities  and  the  pressure  at  the  centre  of  the  cell.  When  the 
continuity  equation  is  added,  the  resulting  5x5  system  can  be  solved  for  each  cell. 

Two  problems  need  attention  with  regard  to  recirculating  gas  flows;  securing  compliance 
with  the  preset  inlet  boundary  conditions  and  selecting  a  discretisation  algorithm  which 
maximises  accuracy  and  secures  convergence. 


3.3.1  The  Momentum  equations  near  inlets 

A  brief  description  of  the  discretised  transport  equations  and  the  method  of  solution  is  given 
below. 

The  discretised  form  of  the  2D  Navier-Stokes  equations  applied  to  the  four  faces  of  a  cell  is 
(from  Refs  [1]  &  [2]): 


APUuiUui+DUvPp  =  ^  {AiUw^i)+ DUwPw -^dpf dr)  DUv+SUUw+SPUv,u^  ,  (S6) 

N,S,E,W 

APU^u^-  DUyPP  ^  53  (AiU^Ui)-DUfPB-\-dpldri\pDV,  +  SUU,  +  SPU,u^,  (87) 

N,S,E,W 

APV.v,  +  DV.pp=  53  {AtV.Vi)  +  DV.ps  +  dp/d^UDV,  +  SUV.  +  SPV,v.,  (88) 

N.S,B,W 

APVpUp  -  WpPP  =  53  {AiV^Vi)~DV,pN  +  dp/d(lDV^  +  SUV^ASPV^v„.  (89) 

N,S,E,W 

(In  the  References  the  pressure  cross-derivatives  dpfdrj  and  dpjd^  are  included  in  the 
source  term) 


As  in  Section  3.2,  DU  and  DV  are  the  cell  face  areas,  and  p  the  static  pressure.  Lower 
case  indices  in''’!. '.ate  cell  faces  and  upper  case  indices  indicate  cell  centres,  except  for  the 
summations  which  are  expanded  as,  for  instance: 


53  (AiUpUi)  ANUpUn  -f  AEDpUe  +  ASUpU,  A  AWUpU^  .  (90) 

N,S,E,W 
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Id  these  equations  APU,  APV  are  the  coefficients  multiplying  u  and  v  as  a  result  of  the 
discretisation  process;  the  subindex  indicates  the  face  at  which  they  are  computed.  Sur¬ 
rounding  velocities  Ui  and  Ui  are  grouped  on  the  RHS,  multiplied  by  similar  coefficients  AiU 
and  AiV,  where  t  stands  for  the  four  surrounding  cells  N,  S,  E  and  W. 

Generalised  coordinates  ((,1;)  are  used,  which  explains  the  pressure  cross-derivatives;  for 
instance  the  first  equation,  which  is  the  f/  or  {  momentum  equation,  contains  a  derivative 
with  respect  to  rf.  Because  boundary-fitted  coordinates  are  often  furly  well  aligned  with 
the  60W,  these  derivatives  are  usually  numerically  small  and  are  not  manipulated  further; 
they  are  later  added  to  the  residuals. 

The  source  terms  are  given  in  a  linear  form  for  the  sake  of  completeness,  comprising  an 
independant  source  term  {SUU,  SUV)  computed  at  the  ceil  face,  and  a  term  linearly 
proportional  to  velocity  {SPU,  SPV).  For  the  Navier-Stokes  equations  SPU  and  SPV 
are  zero. 

The  current  values  of  the  primitive  variables  u*,  «*  and  p*  do  not  satisfy  the  system  of  FDE’s, 
and  a  residual  can  be  computed  equal  to  the  left-hruid  side  minus  the  right-hand  side  of 
each  equation  computed  with  current  values.  Then,  substituting  the  correction  expressions: 


u  =  u*  -1-  u  ,  (91) 

V  =  v'  +  v' ,  (92) 

P  =  P*  +  P  .  (93) 

there  result  the  correction  equations; 

APU^u'^  +  DU^p'p  =  (AiU^i^'i) DU^p'w  +  ^  (94) 

N,S,E.W 

APU, u^-DU,Pp  =  Y,  (AiU^u'i)  -  DU,p'b  +  RUp  ,  (95) 

N,S,E,W 

APV.v,  +  DV.pp  =  Y,  (AiV,v'i)  +  DVps-\-RV.,  (96) 

N,S,E,W 

APV, v^  +  DV,pp  =  (AiV^v^)-DV^p„  +  RV,,  (97) 

N,S,E,W 


where  RU,  RV  are  the  residuals. 


As  mentioned  above,  a  5x5  system  of  ordinary  FDE’s  is  of  ined  (called  the  Block  Implicit 
System)  which  can  be  arranged  as; 

011  0  0  0  ois 

0  (>22  9  0  025 

0  0  033  0  035 

0  0  0  044  045 

.  <>51  052  <*53  <*54  <*S5 
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«-momentum  at  w  , 

t 

F2 

u-momentum  at  p , 

r 

= 

F3 

u-momentum  at  s  ,  (98) 

/ 

U. 

F, 

v-momentum  at  p , 

r 

L  pp  \ 

F, 

Continuity  . 

4r 


•  • 


•  •  • 


•  • 


•  • 


I 


i 

•  I 


Coefficioits  Uij  can  be  identified  by  comparison  with  the  FDE’s;  for  instance,  for  the  u- 
momentum  equation  at  the  to  face; 


On  =  APU^,  (99) 

Ois  =  DU^  ,  (iOO) 

and  for  the  continuity  equation  (  Eq.  (73),  Section  3.2.2),  the  correspondence  is  051  =  oi, 
oss  =  02,  etc. 


This  system  is  solved  for  each  cell.  Then,  all  cells  are  solved  in  one  sweep  of  the  flow  region 
(note  that  the  solutions  are  transported  between  cells  since,  e.g.,  of  one  cell  becomes 
u,„  of  the  adjacent  cell  to  the  right,  and  various  u',  v  and  p  are  placed  on  the  RHS’s  of 
other  cells).  The  whole  field  is  then  swept  repeatedly  until  a  converged  solution  is  found, 
indicated  by  a  reduction  of  the  residuals  to  acceptable  levels. 

New  values  of  the  primitive  variables  are  then  found  as: 


u  = 

u*  +  u'  , 

(101) 

V  = 

u*  +  w' , 

(102) 

p  = 

p'  +  p  . 

(103) 

p  = 

.  .  (P-VH') 

(104) 

This  completes  the  brief  presentation  of  the  discretised  transport  equations  imd  the  method 
of  solution. 

A  particular  problem  euises  near  boundaries,  where  primitive  variables  ue  specified  and  are 
not  to  be  corrected.  For  instance,  if  the  to  face  is  a  solid  wall  the  following  must  be  complied 
with: 


u*  =  “«  =  0  >  (105) 

<  =  ««  =  0.  (106) 

This  is  catered  for  by  making,  in  this  case,  uu  =  1,  uu  =  Fi  =  0.  This  yields  =  0,  and  the 

value  of  Uw  will  be  unperturbed  (i.e.,  will  continue  to  be  zero  as  set  initially).  However,  p'p 

is  now  disconnected  from  the  u-momentum  at  the  to  face,  i.e.,  from  the  prescribed  boundary 
values. 

This  introduces  the  problem  of  prescribed  boundary  values.  The  necessary  and  sufficient 
conditions  at  boundaries  for  compressible  viscous  flows  are  not  theoretically  known,  and 
boundary  conditions  are  usually  specified  in  an  ad  hoc  manner.  Since  the  conservation  of 


momentum  and  energy  equations  are  to  be  satisfied  the  values  of  u,  v,  and  t  are  set,  either 
by  Neumann  (gradient)  or  Dirichlet  (value)  conditions. 

As  the  continuity  equation  must  2dso  be  satisfied  and  u  and  v  have  been  set,  the  last  value  to 
be  specified  is  density  or,  since  continuity  has  been  put  in  pressure-velocity  terms,  pressure. 

For  solid  wall-,  outlet-  and  symmetry-type  boundaries  the  last  condition  is  a  simple  Neu¬ 
mann  condition: 


Q 

—  =  0  (n=  normal  to  boundary) 
On 


(107) 


This  is  effected  by  making  the  pressure  at  the  boundary  equal  to  the  pressure  at  the  first 
inner  cell  centre. 

This,  however,  leaves  the  general  pressure  level  totally  free.  Inlet-type  boundaries  are  the 
only  means  left  to  link  the  pressure  levels  over  the  whole  field  to  the  fluid  inlet  conditions. 
Hence,  a  Neumann  condition  is  applied  in  reverse,  i.e.,  the  boundary  pressure  is  given  and 
the  first  inner  cell  centre  forced  to  this  value.  In  the  correction  system  p  would  be  forced 
to 


P  —  Pioundary  ~  P  •  (1®®) 

This  is  the  usual  procedure,  i.e.,  to  assume  no  streamwise  static  pressure  gradient  at  the 
inlet,  and  becomes  more  satisfactory  as  the  mesh  is  made  finer.  In  the  multigrid  system, 
however,  quite  coarse  meshes  may  be  used,  and  the  assumption  is  not  satisfactory. 

To  achieve  a  better  inlet  boundary  condition  specification,  start  again  from  the  second  law: 


ds  _dt  7  —  1  dp 

Cf~  t  7  P  ' 


(109) 


Since  the  transformation  between  static  and  stagnation  states  is  by  definition  isentropic  this 
may  also  be  written: 


ds  dT  7  -  1  dP 

Cy~  T  'f  P  - 


(110) 


It  is  assumed  that  at  the  inlet  the  flow  will  not  be  reacting,  and  that  the  enthidpy  difference 
between  cells  is  relatively  small,  whence  heat  transfer  by  radiation  and  conduction  between 
cells  at  the  inlet  is  negligible.  Then: 


dT 

T 


=  0. 


(Ill) 


Also,  as  before. 


Substituting  and  manipulating: 


In  generalized  coordinates: 


ds  =  —dr 

pi 


dP  P. 

dr  p 


(112) 


(113) 


dr 

dP 

dx 

dy 


dP  dP 

i 

J\dr)d^  didf!) 
J  dri  dr)  j 


(114) 

(115) 

(116) 

(117) 

(118) 


Substituting  and  manipulating: 


(119) 


If  4  and  T)  are  boundary-fitted  coordinates  and  the  inlet  boundary  has  been  sensibly  chosen 
(i.e.,  with  minimal  divergence  of  the  velocity  vector),  one  of  the  coordinates  (say,  ()  will  be 
substantially  aligned  with  the  flow.  Focusing  then  on  dPId^  (Fig  10): 


dP  P-Pw 
di  ~  Ae/2 


Substituting  and  rearranging: 


Since  time  is  invariant: 


(120) 


(121) 


JA^/U  ~  Ax/u  ~  At  . 


(122) 
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This  yields: 


P  = 


Pw  - 


VdP^ 

2 


1  + 


(123) 


This  equation  links  the  total  pressure  at  the  first  cell  centre  P  to  the  total  pressure  at 
the  inlet  Pw-  There  is  a  total  pressure  drop  due  to  the  work  of  the  viscous  forces  (the 
denominator  is  always  larger  than  unity).  The  second  term  in  the  numerator  represents 
transport  from  adjacent  cells  into  the  control  volume  of  fluid  with  a  different  total  pressure, 
in  the  transverse  direction.  In  the  conventionally  positive  case,  for  example,  with  flow  from 
bottom  to  top,  V/U  >0,  dPfdr)  <  0,  and  there  is  a  net  contribution  which  reduces  the 
total  pressure  loss. 

The  expression  allows  calculation  of  P  and  hence  the  static  pressure,  from  where  the  pressure 
correction  is  calculated  as: 


P  =  PealeulaieJ  —  P"  ■  (124) 

This  is  an  exact  formulation  and  is  hence  applicable  to  the  coarsest  meshes.  In  practice  it  is 
found  sufficient  to  enforce  it  on  the  finest  meshes  and  leave  the  coarsest  meshes  free,  which 
appears  to  accelerate  convergence. 

3.3.2  The  Power-Law  scheme 

This  Section  describes  improvements  made  to  the  procedure  used  to  generate  the  coefficients 
in  the  discretised  transport  equations. 

The  coefficients  of  u  and  v  in  the  transport  FDE’s  are  AW,  AiV  and  APU,  APV,  where 
t  stands  for  the  four  surrounding  centres  E,  W,  N  and  S,  and  the  last  two  coefficients  are 
obtained  by  summing: 


APU  =  (AW), 

(125) 

APV  =  Y,  (^»'')- 

(126) 

ff,E,S,W 


Expressions  for  these  coefficients  vary  according  to  the  differencing  scheme  adopted  to  con¬ 
vert  the  Partial  Differential  Equations  (e.g.,  the  Navier-Stokes  equations)  into  Finite  Dif¬ 
ferences  Equations.  If,  for  instance,  a  Taylor  series  expansion  centred  on  the  cell  face  is 
used,  the  result  is  known  as  the  Central  Differences  scheme,  and  the  coefficients  are,  for 
example 


AEU  =  DX-CX, 

(127) 

AWU  =  DXaCX, 

(128) 

ANU  =  DY-CY, 

(129) 

ASU  =  DY  +  CY, 

(130) 
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where,  disregarding  coordinate  transformation  metrics,  the  Convection  coefficients  CX, 
CY  and  the  Diffusion  coefficients  DX,  DY  are: 


CX  = 

\pUDU, 

(131) 

CY  = 

\pVDV, 

(132) 

DX  = 

VDU 

Ay  ’ 

(133) 

TDV 

DY  = 

Ax 

(134) 

A  detailed  discussion  of  this  scheme  may  be  found  in  [8],  where  it  is  also  shown  that  when 
the  cell  Peclet  number 


m 


•  • 


Pe  = 


2CX 

DX 


exceeds  2  in  absolute  value  the  solution  may  become  unstable. 


(135) 


There  are  a  several  methods  which  may  be  used  to  restrict  the  excursions  in  the  Central 
Difference  scheme  whilst  miuntaining  the  desirable  characteristics  of  the  scheme,  i.e.,  ability 
to  compute  reversing  flows,  and  second-order  accuracy.  Incorporation  of  such  methods  to 
MULTIFL0W-2D  is  not  currently  envisaged  since  the  multigrid  system  allows  the  use  of 
very  fine  grids,  and  gas  turbine  combustion  does  not  involve  relatively  high  gas  velocities. 

Stability  may  be  secured  by  resorting  to  the  Upwind  Differencing  scheme,  where  the 
convection  terms  are  discretized  using  the  principle  that  the  value  of  the  vuiable  at  the 
cell  face  equals  the  value  at  the  upwind  cell  centre,  much  as  was  described  in  Section  3.2.1 
for  density.  Diffusion  terms  are  discretized  as  above.  This  scheme  is  less  accurate,  and 
unsuitable  for  recirculating  flows,  but  may  be  acceptable  since  it  will  be  used  in  the  coarse 
grids,  and  in  regions  of  flow  where  convection  is  dominant. 

In  this  case  the  coefficients  become,  for  example: 


•  • 


•  •  • 


•  • 


•  • 


AEU  =  DX  +  M[-2CX],  (136) 

AWU  =  DX  +  M[-i-2CX]  ,  (137) 

•  • 

where  Af  is  the  MAX  operator  defined  before  (Section  3.2.1).  For  example,  if  CX  >  0 
(conventionally  positive  U),  two  coefficients  CX  are  added  to  the  upwind  side  {AWU)  and 
none  to  the  downwind  side  (AEU),  compared  to  one  each  in  the  Central  Differences 
scheme. 

•  • 
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•  •  • 


The  original  version  of  MULTIFL0W-2D  incorporated  these  considerations  by  means  of 
computing  the  coefficients  in  one  or  the  other  manner  according  to  the  value  of  the  Peclet 
number.  This  is  known  as  the  Hybrid  Scheme  of  discretisation. 

This  was  found  unsatisfactory  in  the  multigrid  environment:  since  each  grid  is  twice  the  size 
of  the  next  finer  grid,  it  is  possible  for  regions  of  flow  to  exist  where  the  cell  Peclet  number 
oscillates  around  2,  and  different  discretization  schemes  are  used  in  contiguous  grids.  In 
these  cases  information  transfers  between  grids  not  only  do  not  help  convergence,  but  may 
lead  to  divergence. 

Patankar  [8]  has  shown  that  an  exact  solution  (at  least  in  one-dimensional  transport)  in¬ 
volves  not  an  on/off  change  between  schemes  but  a  smooth  transition  via  an  exponential 
function  of  the  Peclet  number;  he  has  also  given  an  empirical  approximation  to  the  exact 
solution.  With  this  the  coefficients  can  be  written: 


AEU  =  DX  A{\Pt\)A- M[-2CX\,  (138) 

AWU  =  DX  A{\Pe\)-\-M[-\-2CX],e\.<-,  (139) 


where  the  Power-Law  function  A  is  given  by: 


>l(|Pc|)  =  A/((l-i-giy].  (140) 

This  has  been  incorporated  into  the  code. 

The  Power- Law  has  also  been  applied  to  the  generalised  transport  equations  for  scalars 
such  as  enthalpy  and  mixture  fraction,  and  particularly  to  the  K  —  t  turbulence  model.  The 
modification  requires  the  inclusion  of  Prandtl  numbers  a  so  that  the  diffusion  coefficients 
become  DXjff,  DY/a,  and 


Pe  = 


2<tCX 

DX 


etc. 


(141) 


4.  APPLICATION  EXAMPLES 
4.1  Jet  in  a  sudden  expansion 


The  first  example  is  chosen  to  show  the  ability  to  compute  strongly  recirculating  flows.  A 
circular  jet  0.12  m  in  diameter  enters  a  circular  duct  0.6  m  in  diameter  and  2.4  m  long. 
Inlet  conditions  are; 


•  Reference  inlet  velocity  Uq=  200  m/sec, 

•  v=  0, 
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•  p=  104429  Pa, 

•  r=  300  K. 


A  l/7th-power  velocity  distribution  law  is  assumed  at  the  jet  inlet,  formulated  to  yield  the 
same  volumetric  flow  rate  as  Uq.  Hence,  inlet  axial  velocity  at  a  radius  r  is  given  by: 

„  120/.  r\‘/^ . .  „„„ 


120  / 

^  =  with  0.06. 


Inlet  turbulence  is  set  at  5%  of  U  for  all  components: 


^  =  3^0.05//)^. 


Hence  static  temperature: 


Whence  density: 


t  =  T  -  ^  ,  with  Cp=  1000  mVsec^K 


p  =  p/Rt  ,  with  R=  286.7  m^sec^K  . 


Turbulence  scale  length  is  set  at  L  =  R/2.  Hence,  with  Cd=  0.09 


With  the  above,  inlet  conditions  at  the  centreline  2U’e: 


U  —  245m/sec , 

K  =  225m*/sec* , 
t  =  269.8K, 
p  =  1.3485kg/m® , 

£  =  18485m*/sec^  , 

M  =  0.75. 

The  finest  computational  mesh  has  240x32  cells,  and  there  are  three  lower  grid  levels,  i.e., 
coarse  grids  of  30x4,  60x8  and  120x16  cells. 

Computation  results  are  shown  in  Figures  11  to  14.  Fig  11  shows  the  upper  half  of  the  duct 
with  flow  trajectories  obtained  by  integration  of  the  velocity  vectors.  The  centre  jet  enters 
at  left  with  a  0.06  m  radius. 


•  • 


Fig  12  is  a  close-up  view  of  the  region  near  the  inlet,  where  a  small  counter-rotating  vortex 
is  visible  at  the  corner,  to  give  an  indication  of  the  fineness  of  detail  obtainable. 

Fig  13  shows  contours  of  turbulent  kinetic  energy  x;  note  the  closely  packed  contour  lines 
between  the  centreline  and  A  at  x  ~  0.5.  Turbulence  intensities  can  be  computed  from  k 
as: 

(147) 

Expressed  as  a  fraction  of  centreline  velocities,  turbulence  intensities  range  from  13%  near 
the  inlet  to  71%  near  the  outlet. 

Turbulent  kinetic  energy  is  also  plotted  as  a  3D  carpet  plot  in  Fig  14,  with  the  jet  entering 
from  right  to  left  at  the  origin  (0,0).  This  representation  evidences  the  horseshoe  pattern 
of  the  peak  TKE,  which  is  not  obvious  in  Fig  13,  and  was  only  suggested  by  the  packed 
contour  lines. 

The  code  performance  is  deemed  satisfactory,  in  that  it  yields  plausible  flow  structures  with 
the  expected  features  (such  as  recirculation  zones  and  shear  layers),  and  fine  detail. 

4.2  Axisymmetric  Contraction 

The  second  example  is  a  design  verification  exercise  of  a  contracting  circular  duct.  Duct 
diameters  are  0.525  m  at  the  inlet  and  0.135  m  at  the  outlet;  duct  length  is  0.62  m.  The 
duct  profile  was  designed  by  others,  and  it  is  desired  to  verify  that  this  shape  can  accelerate 
a  flow  of  air  by  approximately  an  order  of  magnitude  without  risking  separation  and  with 
uniform  discharge  velocity. 

Inlet  conditions  are  iteratively  adjusted  to  obtain  M  ca  0.8  at  the  outlet  and  a  properly 
expanded  flow  (outlet  static  pressure  approximately  atmospheric).  Inlet  velocity  is  a  nominal 
13  m/sec,  with  a  l/7th-power  radial  distribution  and  5%  turbulence.  Three  mesh  levels  were 
used  with  17x4,  34x8  and  68x13  cells. 

Computation  results  are  shown  in  Figures  15  to  19.  Fig  15  shows  true  streamlines  (surfaces 
enclosing  a  given  mass  flow,  as  distinct  from  the  flow  trajectories  of  Fig  11).  Thickening 
of  the  90%-100%  layer  in  the  concave  region  is  modest,  indicating  a  low  probability  of 
separation.  There  are  no  indications  of  reversed  or  separated  flow. 

Contours  of  axial  velocity  (Fig  16)  are  nearly  circular  in  the  approach  region,  which  is  taken 
as  an  indication  of  a  satisfactory  design  in  that  fluid  parcels  at  similar  distances  from  the 
outlet  possess  the  same  velocity.  The  exit  velocity  distribution  is  quite  uniform  (Fig  17), 
and  the  exit  Mach  number  is  approximately  0.8. 

Turbulence  intensity  is  low  throughout:  from  the  inlet  value  of  5%  it  drops  to  approximately 
2%  at  the  outlet  (Fig  18). 

Fig  19  shows  that  the  radial  static  pressure  gradient  at  the  outlet  is  low,  indicating  a 
substantially  parallel  outlet  flow. 

Detuled  analysis  of  these  numerical  simulation  results  and  comparisons  with  the  physical 
experiment  will  be  reported  elsewhere.  Nevertheless  it  may  be  said  that  the  simulation 
shows  the  duct  design  to  be  satisfactory. 
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5.  CONCLUSIONS 


A  computer  program  named  MULTIFL0W-2D  has  been  acquired  for  use  at  ARL.  This 
program,  which  is  capable  of  computing  numerical  simulations  of  internal  reactive  flows, 
has  been  thoroughly  analysed  and  reflned.  Important  improvements  relevant  to  gas  turbine 
combustion  have  been  developed  and  incorporated,  such  as  the  ability  to  compute  flows 
with  variable  density  and  also  non-isentropic  flows.  Code  development  has  reached  the 
stage  where  it  can  be  applied  with  confidence  to  practical  problems. 

Until  now  development  has  been  focused  on  the  aerothermodynamics  aspects  of  the  code. 
Development  of  the  reactive  flow  facilities  (chemistry,  kinetics)  should  now  proceed  as  is 
planned.  As  a  matter  of  detul,  it  is  intended  to  reassess  the  possible  incorporation  of 
limiting  schemes  during  later  development  stages,  to  dispense  with  the  Power-Law  scheme. 
This  may  be  appropriate,  for  instance,  if  a  need  arises  to  use  higher  order  discretisation 
algorithms. 
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Axisymmetric  expanding  duct  with  2x1  coarse  grid 
and  next  two  fine  grids  (original  scheme) 

FIG  1 


Axisynunetric  expanding  duct  with  8x4  fine  grid 
and  next  two  coarse  grids  (revised  scheme) 
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Upwind- biased  second-order  interpolation 
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Axisymmetric  contraction 
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A  computational  flmd  dynamics  code  named  MULTIFLOW-2D  has  been  acquired  and  is  being  put  into 
service  for  the  numerical  simulation  of  gas  turbine-type  combustors. 

Certain  areas  of  the  code  were  found  to  be  in  need  of  improvement  for  this  application;  these  are: 

•  The  grid  generation  algorithm 

•  The  manner  of  discretisation  of  the  continuity  and  momentum  equations,  and 

•  The  treatment  of  the  inlet  boundary  conditions. 
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